%% Table C1: Labor Market Skills: Moments and Parameter Estimates

%% Code for Solving the Parameters of a Two-Skills Lognormal Roy model for Brazil
%   --- Inferring the means, stds and covariances of Skills
clear all
close all
warning off

%% Number of schooling groups; it can be NE=4, NE=5 (or NE=6?)
        NE = 5;  

%% Earnings data for Brazil, 2000 (from Income_Education_Skills Excel file)        
    %  columns:   (mean, stds of earnings) (mean,std)_r, (mean,std)_n ; 
    %  rows:      education groups


%% NE=5    
% Mean (columns 1 & 3) and standard deviation (columns 2 & 4) of earnings
% by education levels (lines) and occupation (routine: columns 1&2;
% cognitive: columns 3&4).
MU_SIGMA_E0 = [ 6.216	0.666	6.652	1.024
                6.672	0.726	7.415	1.081
                7.007	0.776	7.664	1.016
                7.338	0.821	7.429	1.087
                7.82	0.824	8.295	0.955];
         
% fraction of workers, columns: r,n; rows: e          
F_E        =  [ 0.993774689	0.006225311
                0.976544936	0.023455064
                0.958710794	0.041289206
                0.844586587	0.155413413
                0.520274578	0.479725422  ];

%% Normalization, so that the mean for earnings of e=1 in r is 1;
        MU_SIGMA_E      = MU_SIGMA_E0;
        MU_SIGMA_E(:,1) = MU_SIGMA_E0(:,1)-MU_SIGMA_E0(1,1)-1/2*MU_SIGMA_E0(1,2)^2;
        MU_SIGMA_E(:,3) = MU_SIGMA_E0(:,3)-MU_SIGMA_E0(1,1)-1/2*MU_SIGMA_E0(1,2)^2;
    
            
%% Inputs (Observations): Fraction, and means and of stds of log-earnings (5 objects)
for e=1:NE
    fr        =    F_E(e,1);
    m_y_r     =    MU_SIGMA_E(e,1);
    m_y_n     =    MU_SIGMA_E(e,3);
    s_y_r     =    MU_SIGMA_E(e,2);
    s_y_n     =    MU_SIGMA_E(e,4);
    OBS = [fr, m_y_r, m_y_n, s_y_r, s_y_n];
    
%% Initial Estimates: means, variances and covariance of skills (5 objects)
    mu_r    =  m_y_r;
    mu_n    =  m_y_n;
    S_rr    =  s_y_r;
    S_nn    =  s_y_n;
    S_rn    =   0;
    PARAMS0 = [mu_r, mu_n, S_rr, S_nn, S_rn];
  
%% using fsolve to get an initial condition for fminsearch
  
  fun1=@(params) RoyParams(params,OBS);
            options = optimset('TolFun',1e-10);
            options = optimset('TolX', 1e-10);
            options = optimset('MaxIter',50000);
            options = optimset('MaxFunEvals',50000);
           [log_norm_params1(e,:),fval_solve(e),exitflag_fsolve(e),output]  = fsolve(fun1,PARAMS0,options);
           PARAMS0= log_norm_params1(e,:);  % using fsolve to get an initial condition for fminsearch
           [log_norm_params(e,:), fval_search(e),exitflag_fmin(e),output]   = fminsearch(fun1,PARAMS0,options);
        


end

disp([log_norm_params])

for e=1:NE,
    mu_r_e(e,1)=log_norm_params(e,1);
    mu_n_e(e,1)=log_norm_params(e,2);
    sigma_r_e(e,1)=log_norm_params(e,3);
    sigma_n_e(e,1)=log_norm_params(e,4);  
    rho_rn_e(e,1)=log_norm_params(e,5)/( log_norm_params(e,3)*log_norm_params(e,4) ); %^.5;
end

%% Normalize so that the mean of the skills of lowest e in r is 1
Moments_Skills_e0           = [mu_r_e,mu_n_e,sigma_r_e,sigma_n_e,rho_rn_e];
Moments_Skills_e            =  Moments_Skills_e0;
Moments_Skills_e(:,1:2)     =  Moments_Skills_e0(:,1:2)-Moments_Skills_e0(1,1)-1/2*Moments_Skills_e0(1,3)^2;


disp('Covergence? exitflag_fmin=1')
disp(exitflag_fmin)
disp('[mu_r,mu_n,sigma_r,sigma_n,rho_rn]')
disp(Moments_Skills_e)


